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ABSTRACT 

A numerical scheme that incorporates a thermal leakage injection model into a combined gas dynamics and 
cosmic ray (CR, hereafter) diffusion-convection code has been developed. The hydro/CR code can follow in a 
very cost-effective way the evolution of CR modified planar quasi-parallel shocks by adopting subzone shock- 
tracking and multi-level adaptive mesh refinement techniques. An additional conservative quantity, S = P g / p^ s ~ l , 
is introduced to follow the adiabatic compression accurately in the precursor region, especially in front of strong, 
highly modified shocks. The "thermal leakage" injection model is based on the nonlinear interactions of the 
suprathermal particles with self-generated MHD waves in quasi-parallel shocks. The particle injection is followed 
numerically by filtering the diffusive flux of suprathermal particles across the shock to the upstream region ac- 
cording to a velocity-dependent transparency function that controls the fraction of leaking particles. This function 
is determined by a single parameter, e, which should depend on the strength of postshock wave turbulence, but 
is modeled as a constant parameter in our simulations. We have studied CR injection and acceleration efficien- 
cies during the evolution of CR modified planar shocks for a wide range of initial shock Mach numbers, Mo, 
assuming a Bohm-like diffusion coefficient. For expected values of e the injection process is very efficient when 
the subshock is strong, leading to fast and significant modification of the shock structure. As the CR pressure 
increases, the subshock weakens and the injection rate decreases accordingly, so that the subshock does not dis- 
appear. Although some fraction of the particles injected early in the evolution continue to be accelerated to ever 
higher energies, the postshock CR pressure reaches an approximate time-asymptotic value due to a balance be- 
tween fresh injection/acceleration and advection/diffusion of the CR particles away from the shock. In the strong 
shock limit of Mo > 30, the injection and acceleration processes are largely independent of the initial shock Mach 
number for a given e, while they are sensitively dependent on Mo for Mo < 30. We conclude that the injection 
rates in strong parallel shocks are sufficient to lead to rapid nonlinear modifications to the shock structures and 
that self-consistent injection and time-dependent simulations are crucial to understanding the non-linear evolution 
of CR modified shocks. 

Subject headings: acceleration of particles-cosmic rays- hydrodynamics- methods:numerical 



1. INTRODUCTION 

Two decades ago it was recognized that CRs are proba- 
bly produced very efficiently via diffusive shock acceleration 
(DSA) in ubiquitous astrophysical shocks, (for early reviews 
see, e.g., Drury 1983; Blandford and Eichler 1987; Berezhko 
and Krymskii 1988). After the initial successes of the simple 
concept that the particles can gain energy while temporarily 
trapped in the converging flows around a shock, it was quickly 
realized that the full DSA treatment requires one to consider 
the complex nonlinear interactions between energetic particles, 
resonantly scattering waves and the underlying plasma (Malkov 
and Drury 2001). One of the important aspects of those interac- 
tions is the injection of suprathermal particles into the CR pop- 
ulation at shocks. According to quasi-linear theory as well as 
plasma simulations of strong quasi-parallel shocks, the stream- 
ing motion of superthermal particles against the background 
plasma can induce wave generation leading to strong down- 
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stream MHD waves that scatter particles and inhibit the parti- 
cles from leaking upstream (e.g., Bell 1978; Quest 1988). Par- 
ticles in or close to the thermal population are especially re- 
stricted in this way. As a consequence only a small fraction 
of suprathermal particles can swim upstream against the wave- 
particle interactions in the plasma flow and be injected into the 
higher energy CR population to be further accelerated via the 
Fermi process. The injection process and its efficiency control 
the amplitude of the CR population and hence the degree of 
shock modification. 

Recently, significant progress in understanding this injec- 
tion process in parallel shocks has been made through self- 
consistent, analytic, nonlinear calculations by Malkov and Volk 
(1995) and Malkov (1998). The resulting theory has only one 
parameter; namely, the intensity of the downstream waves, and 
that is tightly restricted, both by the theory and by comparison 
with hybrid plasma simulations. By adopting Malkov's ana- 
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lytic solution, we have developed a numerical treatment of this 
injection model and incorporated it into a combined gas dynam- 
ics and CR diffusion-convection code (Gieseler et al. 2000). 
According to the Gieseler et al. (2000) simulations, the injec- 
tion process seemed to be self -regulated in such a way that the 
injection rate reaches and stays at a nearly stable value after 
quick initial adjustment, but well before the CR shock reaches a 
steady state structure. Gieseler et al. (2000) found about 10~ 3 of 
incoming thermal particles to be injected into the CRs, roughly 
independent of Mach numbers. However, due to severe compu- 
tational requirements associated with the need to resolve struc- 
tures down to the physical shock thickness, those simulations 
were carried out only until the characteristic maximum momen- 
tum of (p max /m p c) ~ 1 was achieved. Since strong shocks were 
still evolving at the end of the simulations, the time-asymptotic 
limit could not be estimated for either the CR acceleration effi- 
ciency or the CR spectrum. 

Unlike ordinary gas shocks, the CR shock is a collisionless 
structure, and includes a wide range of length scales associated 
not only with the dissipation into "thermal plasma", but also 
with the nonthermal particle diffusion process. Those are char- 
acterized by the so-called diffusion lengths, DdiffQ?) = k(p)/u, 
where n(p) is the spatial diffusion coefficient for CRs of mo- 
mentum p, and u is the characteristic flow velocity (Kang 
and Jones 1991). For strong scattering of suprathermal parti- 
cles Ddiff may not greatly exceed the physical dissipative, or 
"gas" shock thickness . Accurate solutions to the CR diffusion- 
convection equation require a computational grid spacing sig- 
nificantly smaller than D^ff, typically, Ax <~ O.lDdifK/?). On 
the other hand, for a realistic diffusion transport model with a 
steeply momentum-dependent diffusion coefficient, the highest 
energy, relativistic particles have diffusion lengths many orders 
of magnitude greater than those of the lowest energy particles. 

To follow the acceleration of highly relativistic CRs from 
suprathermal energies, all those scales need to be resolved nu- 
merically. However, the diffusion and acceleration of the low 
energy particles are important only close to the shock owing 
to their small diffusion lengths. Elsewhere, they are effectively 
advected along with the underlying gas flow. At higher ener- 
gies the needed resolution is less severe, and on scales larger 
than ~ £>diff(/?max) the bulk flow and the nonthermal particles 
decouple, so resolution requirements are controlled by what- 
ever factors are necessary to define the broader flow properties. 

Thus it is necessary to resolve numerically the diffusion 
length of the particles only around the shock. So, in Kang et 
al. (2001) we first implemented a shock tracking scheme (LeV- 
eque and Shyue 1995) to locate the shock position exactly and 
then refine the grid resolution only around the shock by apply- 
ing multi-levels of refined grids (Berger and LeVeque 1998). 
The main properties of this code are: 1) the shock is tracked as 
an exact discontinuity, 2) a small region around the shock is re- 
fined with multi-level grids, 3) it is very cost-effective in terms 
of computational memory and time. 

In the present contribution, we have studied the CR injec- 
tion and acceleration during the evolution of modified planar 
shocks by implementing the numerical method for the thermal 
leakage injection model of Gieseler et al. (2000) into an en- 
hanced version of the CR/AMR hydrodynamics code, which 
we name CRASH (Cosmic Ray Amr SHock) code. With our 
new CRASH code we were able to calculate the CR injection 
and acceleration efficiencies with a Bohm-like diffusion coef- 
ficient for higher particle energies (p/m p c » 1) in shocks over 



a wide range of initial Mach numbers. In the following section 
the basic equations solved in the simulations are presented. We 
describe the thermal leakage injection model and its numeri- 
cal implementation in our CR transport code in §3. We then 
outline our numerical methods in the CR/AMR hydrodynamics 
code in §4. In §5 we present and discuss our simulation results, 
followed by a summary in §6. Finally we discuss some test 
calculations in the Appendix. 

2. BASIC EQUATIONS 

We solve the standard gasdynamic equations with CR pres- 
sure terms added in the conservative, Eulerian formulation for 
one dimensional plane-parallel geometry: 



dp | d(up) = 
dt dx 

d(pu) + d(pu 2 + P g + P c ) 



dt 



dx 



■0, 



d(pe g ) d(pe g u + P s u + P c u) 



= -L(x,t), 



(1) 



(2) 



(3) 



dt dx 
where P g and P c are the gas and the CR pressure, respectively, 
e g = Pg/ P(lg~ i)+u 2 /2 is the total energy density of the gas per 
unit mass and the rest of the variables have their usual mean- 
ings. The injection energy loss term, L(x,t), accounts for the 
energy of the suprathermal particles injected to the CR compo- 
nent at the subshock. L is nonzero only at the subshock tran- 
sition, as explained at the end of §3.2 and in §4.1. CR inertia 
is neglected, as usual, in such computations. A basic outline of 
the computational hydrodynamical scheme and its adaptations 
for this problem are given in §4. 1 . 

We do include here one unusual adaptation in the CRASH 
code that significantly improves treatment of the CR precursor 
region. In conventional conservative, gas dynamical numerical 
methods, the thermal energy is calculated by subtracting the ki- 
netic energy from the total gas energy e g . For highly supersonic 
flows where the kinetic energy is much larger than the thermal 
energy, however, the numerical errors in calculating the kinetic 
gas energy could be much larger than the thermal energy it- 
self. Especially in strong CR modified shocks where the gas 
flow is compressed through the precursor but remains highly 
supersonic, it is difficult to follow this preshock compression 
accurately. In order to alleviate this difficulty, following the 
approach first introduced by Ryu et al. (1993), we adopt an ad- 
ditional conservative quantity, the "Modified Entropy", 



s = p s /p 



7«- 



(4) 



related to the gas entropy per unit volume. The modified en- 
tropy satisfies the following conservation equation (Ryu et al. 
1993): 

which is valid outside the dissipative subshock; i.e., where the 
gas entropy is conserved. Hence the modified entropy equation 
is solved outside the subshock, while the total energy equation 
is applied across the subshock. As we will discuss below, the 
injection process depends on the properties of the subshock and 
so on the precursor compression. Our new hydro code with 
the modified entropy equation can accurately follow the adia- 
batic compression through the precursor and provides a better 
method to follow the injection process. 
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The diffusion-convection equation for the pitch angle aver- 
aged CR distribution function f(p,x,t) (e.g., Skilling 1975) is 
given by 

df df i du df a df 

af ax 3 <xc op Ox ox 
and k(x,p) is the diffusion coefficient. For convenience we al- 
ways express the particle momentum, p in units m p c. As in our 
previous studies, the function g(p) = p A f(p) is evolved instead 
of f(p) and y = ln(p) is used instead of the momentum variable, 
p for that step. This leads the equation solved to take the form 
(Rang and Jones 1991; Gieseler et al. 2000): 



dg dg 1 du dg A .d dg 
Ot ox 3 ox Oy ox ox 



(7) 



Our numerical approach to this equation is outlined in §4. 1 with 
additional details provided in the cited literature. 

3. THERMAL LEAKAGE INJECTION MODEL 

3.1. Transparency Function 

In the "thermal leakage" injection model, most of the down- 
stream thermal protons are confined by nonlinear waves and 
only particles well into the tail of the Maxwellian distribution 
can leak upstream across the subshock. Since any escaping par- 
ticles have to swim against the advection flow downstream, the 
breadth of the thermal velocity distribution relative the down- 
stream flow velocity in the subshock rest-frame is central to the 
injection problem. In order to model this injection process in 
Gieseler et al. (2000) we adopted a "transparency function", 
r esc , which expresses the probability that supra-thermal parti- 
cles at a given velocity can leak upstream through the mag- 
netic waves, based on non-linear particle interactions with self- 
generated waves (Malkov 1998). In this scheme, the trans- 
parency function is approximated by the following functional 
form, 

r eS c(e, v/u d ) = H[i,-(l + e)](l-^y l (l-^J 

•exp{-[{3-(l + e)r 2 } , (8) 

which depends on the ratio of particle velocity, v, to down- 
stream flow velocity in the subshock rest-frame, uj = u s /r su b, 
where u s and r su b are the subshock velocity with respect to the 
immediate upstream flow and the subshock compression factor, 
respectively. Here H is the Heaviside step function. The param- 
eter, e = Bq/B±, is defined in Malkov (1998), and measures the 
ratio of the amplitude of the postshock MHD wave turbulence 
B± to the general magnetic field aligned with the shock normal, 
Bo- Here v = e v/ud is the normalized particle velocity. 

The lower right panel of Fig. 1 shows the Maxwellian distri- 
bution function, gu(p) = /m(p)p 4 for the subshock Mach num- 
ber M s = 100, and T esc (p) as a function ofv/uj for e = 0.3 (solid 
line) and e = 0.2 (dashed line). The function g M is plotted for 
P = u s /c = 0.01. The function T esc plotted as a function of v/ud 
is independent of (3. The plot for T esc (/?) shows that the injec- 
tion occurs predominantly for the particles with 5ud < v < 20ud 
where < r esc (/?) < 1 . The lower velocity particles with r esc ~ 
represent the gas particles, while the higher velocity particles 
with T esc Ri 1 represent the CR particles. So the particles are 
injected over a range of particle momenta rather than a singly 
chosen injection momentum. Thus the transparency function 
eliminates the need to set up the lowest momentum boundary 
that arbitrarily and abruptly distinguishes thermal plasma and 



accelerated particles. If one takes a step function (i.e., T esc = 
for p < /?; n j and r esc = 1 for p > p^) rather than a smooth func- 
tion for the transparency function, the injection model becomes 
the same as the earlier injection models where the particles 
are injected at an injection momentum, p; n j, which is then the 
momentum boundary between thermal and CR particles (e.g., 
Berezhko et al. 1995; Rang and Jones 1995). 

The only free parameter of the adopted transparency function 
is e and it is rather well constrained, since 0.3 < e < 0.4 is in- 
dicated for strong shocks (Malkov and Volk 1998). Due to the 
exponential cut off in a thermal velocity distribution, however, 
the injection rate depends sensitively on the value of e, since the 
location in the Maxwellian tail distribution where dr esc / dp=/0 
depends on e (see the lower right panel of Fig. 1). Since the 
wave generation is relatively less efficient for weak shocks, re- 
sulting in larger values of e, this parameter should be depen- 
dent on the subshock Mach number. However, since it is not 
certain precisely how e should vary as the subshock becomes 
weaker with time, we will keep the value of e constant within 
each simulation except for a few models where e is changed 
during the simulation to avoid a numerical problem during ini- 
tialization. We expect the injection rate would be higher if we 
let e increase as the subshock weakens, compared with the cal- 
culation with a constant value. Although larger values of e nat- 
urally lead to higher injection rate and higher acceleration ef- 
ficiency, as we will show in $ 5.1.2, the net enhancement of 
CR energy would be minimal (a few %), since the CRs injected 
earlier are dominant over the fresh particles injected recently at 
much weaker shocks with a much smaller injection rate. An- 
other degree of freedom in this model is the detailed functional 
form of the transparency function. The form given by equation 
(8), which is an approximation of the representation given in 
Malkov and Volk (1998), behaves like a smoothed step func- 
tion which selects the range of injection momenta for a given 
value of e. With a modification in the shape of the transparency 
function (due to a more accurate approximation), we expect the 
resulting injection rate may not change significantly, as long as 
the adopted T esc has the characteristics of a smoothed step func- 
tion and picks the similar momentum range for the injection 
pool. 

3.2. Numerical Scheme for Injection Model 

In Gieseler et al. (2000) we proposed a new numerical 
scheme for injection to be used in kinetic equation simula- 
tions of diffusive shock acceleration. The scheme emulates 
thermal leakage injection through the transparency function 
given by equation (8) for particle leakage across a shock dis- 
continuity. The total particle distribution function consists of 
Maxwellian and CR components (g = g M + gctd, and there is 
no discrete momentum boundary between the two components. 
In our scheme the Maxwellian distribution, gw(p,x,t) is cal- 
culated from the density and temperature of the local gas at 
each time step from gas dynamical data, since it cannot be fol- 
lowed by the diffusion-convection equation (6). The suprather- 
mal particles in this "injection pool" should have a significant 
anisotropic component, so the distribution there cannot be fol- 
lowed by the diffusion-convection equation either. In point of 
fact, one cannot cleanly demark particles as thermalized and 
fluid-like or nonthermal and diffusive, yet economical compu- 
tational schemes exist only in those two limits. That is the very 
reason we need an injection model as represented here through 
the transparency function. 
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We briefly explain how the distribution function is integrated 
in a way that can implement the thermal leakage process in our 
CR transport code. 

1) The total particle distribution is updated in distinct "advec- 
tion and diffusion" steps for the momentum bins with T esc > 0: 
g" — > g n+ \ where g n is the value at the previous time step, and 
g" +1 is an intermediate value for g n+x . For thermal particles with 
Tesc = 0, g n+1 = #m 1 • Then the CR distribution up to this inter- 
mediate step is given by g^ = g n+l -g^ 1 . 

2) The net increase of CR particles, (g^ ~<?cr)' calculated 
from the diffusion-convection equation (6) should be reduced 
according to the transparency function in the upstream region. 
So an additional "leakage" step is applied to the CR distribution 
function in the upstream region as 

Scr=Scr + w(0-[scr-Scr]- (9) 
In other words, we first estimate the number of suprather- 
mal particles that cross the shock according to the diffusion- 
convection equation, and then we allow only some fraction of 
the combined advective and diffusive fluxes to leak upstream 
with the probability prescribed by T esc . Obviously the CR par- 
ticles with T esc = 1 are not affected by the leakage step. 

We then need to estimate the gas energy loss rate at the 
subshock, L(x,t), due to particle leakage in order to conserve 
total energy. Although the injection process is realized nu- 
merically via the leakage of suprathermal particles to the up- 
stream region (i.e., diffusion in space), we formulate the en- 
ergy loss rate due to injection in terms of the particle accel- 
eration in momentum space, because leaking particles experi- 
ence an adiabatic energy increase upon crossing the subshock. 
The particle number flux crossing the momentum boundary 
at p due to the adiabatic change of the momentum through 
the flow compression/expansion is given by r esc (p)Q(p), where 
Q(p) = -(4ir /3)p 3 f(p)(du /dx), since T esc (p) gives the probabil- 
ity for the suprathermal particles to cross the shock. Then the 
net change of the particle number density per unit time in the 
momentum bin bounded by p and p + dp can be written as 

dN p = T esc (p + dp)Q(p + dp)-T esc (p)Q(p) 



-z—Q(p)dp + T esc (p)—dp 
op up 



(10) 



If the transparency function is a step function (i.e., dr esc /dp = 
d(Pmi) X the first term of equation (10) gives the injection rate 
at p = pi n j, while the second term represents the acceleration 
of CRs for p > pi n j. For a smooth transparency function, the 
second term contains the acceleration of both CRs and the 
suprathermal particles, since the momentum boundary between 
thermal and CR population is soft and f(p) = /qr + fiw for the 
momenta in the injection pool (i.e., < T esc (p) < 1). The in- 
jection rate due to the second term in equation (10) that corre- 
sponds to the acceleration of suprathermal particles should be 
small because of exponential decrease of /m in this momentum 
range. As shown in Fig. 1, dr esc /dp, for the transparency func- 
tion adopted here behaves almost like the delta function, and 
so the first term of equation (10) dominates for the momenta in 
the injection pool. So we take only the first term to estimate 
approximately the energy loss rate according to 



2tt 



du 



L(x,0 w -— m p c — 
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dr esc (p,t) 5 

P f(P,x,t) dp, 



dp 



(11) 



where we have used the nonrelativistic approximation for the 
kinetic energy of injected particles, jtn p c 2 p 2 . We remind the 
readers that p is expressed in units of m p c. The term L(x,t)At 



is then subtracted from the gas thermal energy in the immediate 
postshock region, and only there. No matching term is nec- 
essary in updating the CRs, since L is introduced explicitly to 
account for evolution of f(x,p,t) in the "injection pool" where 
/ is in transition between thermal and nonthermal forms. Since 
the fraction of injected particles is very small, the energy ex- 
changed through L is, in any case, small, so that it has little 
direct influence on the flow. 

Since the diffusive flux of particles depends on the particle 
diffusion length, the CR injection rate numerically realized in 
our scheme depends on the ratio of Auff(/>inj)/ Ax at the finest 
grid. However, the numerical results seem to be converged at a 
spatial resolution of Ax ~ Auff(Pinj)- 

3.3. Dependence of Injection Rate on the Subshock Mach 

Number 

As we will show below, the ensuing injection rate based on 
our thermal leakage model depends largely on the subshock 
Mach number, M s , and the inverse wave-amplitude parameter, 
e, because the transparency filter is a function of v = e v/ud- 
Here we will explore the consequences of this dependence. The 
thermal peak velocity of the Maxwellian distribution of the im- 
mediate postshock gas relative to the mean flow is defined as 
t>th = 2yJkftT2/m v , where T2 is the postshock temperature, kg is 
the Boltzmann constant, and m p is the proton mass. The mean 
molecular weight is assumed to be one in this consideration and 
also in our numerical calculations. The ratio of v^/ud is then 
determined only by M s as follows: 



Ud 



'\2(5M 2 -\) 
T (M? + 3) ' 



(12) 



where M s is the subshock Mach number and we assume j g = 
5/3. Since the injected particle flux is approximately propor- 
tional to (dr esc / dp) f(p)p 3 dp (see equation [11]), we define the 
"injection velocity", v\ n j as 



J v (dT esc /dp)f M (p)p 3 dp 


mj _ 00 ' 

/ d(T etiC /dp)f M (p)p 3 dp 




Vini = 



(13) 



where /m is the Maxwellian velocity distribution of the im- 
mediate postshock gas. One can easily show that the ratio of 
Vwj/ud is a function of M s and e only (see Fig. 1). Since in 
the model used here the rate of particle injection into the Fermi 
process is proportional to the convolution of dr esc /dp with the 
Maxwellian tail distribution, the injection rate is determined 
mainly by the ratio of finj/fth. rather than the velocities them- 
selves. This ratio depends on only M s and e. We plot the values 
of finj/fth m Fig- 1 for e = 0.2-0.35. For strong shocks M s > 10 
this ratio becomes independent of M s , while it increases for 
smaller M s . To illustrate, for e = 0.3; Wi n j/w th = 1 .87 forM, > 10, 
while Winj/wth — * 2.4 as M s — > 2. Thus, the injection rate be- 
comes independent of the subshock Mach number for strong 
shocks. Larger values of v^/v t \, for smaller M s means that the 
injection process is less efficient for weak shocks. This ratio 
is also larger for smaller e, since the leakage is prohibited by 
stronger wave fields. Thus, injection is more strongly inhibited 
by smaller e (stronger waves). 

Since the injected particles are nonrelativistic, the ratio 
Winj/wth, is in fact the same as the parameter c\ = p\ajpth de- 
fined as a constant injection parameter in an earlier version of a 
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thermal leakage model by Kang and Jones (1995) where values 
in the similar ranges (i.e., c\ = 1.9-1.95) were adopted. 

As shown above, in a "thermal leakage" injection model 
where the leakage probability is determined by the ratio of the 
particle velocity to the advection velocity of the downstream 
flow relative to the subshock, we can expect the following : a) 
the injection rate becomes independent of the subshock Mach 
number forM 5 J> 10. b) the injection is less efficient for weaker 
shocks (M s < 10) compared to stronger shocks, c) in the time 
evolution of a CR modified shock the injection rate decreases 
as the subshock weakens by way of precursor compression. 

4. NUMERICAL METHODS AND MODEL PARAMETERS 

4. 1 . CRASH: A CR/AMR Hydrodynamics Code with Modified 

Entropy 

Equations (l)-(3) and (5) are solved in two steps on an Eu- 
lerian grid: a gasdynamic step without CR terms followed by 
a CR-dynamic step, as outlined below. Then equation (7) is 
solved in two steps to update the CR distribution and allow 
for injection at shocks. The advection term of equation (7) is 
solved by the wave-propagation method, as for the gasdynamic 
variables, except that only the entropy wave applies. Then the 
diffusion term is solved by the semi-implicit Crank-Nicholson 
scheme, as described in Kang and Jones (1991). 

A full description of the combined gas dynamic/CR algo- 
rithm without equation (5) can be found in Kang et al. (2001). 
We refer readers to that work for most details and a discussion 
of various algorithm tests. Since we introduce two new features 
in the present work; namely an improved treatment of the shock 
precursor and our new "thermal leakage injection" model, we 
briefly outline the basic scheme here. Some further, relevant 
code tests are described briefly in the Appendix. 

In the gas dynamic step, the hydrodynamic conservation 
equations are solved by the explicit, LeVeque-Shyue Riemann- 
solver-based wave-propagation algorithm, temporarily ignor- 
ing the CR terms in eq. 2 and 3. A nonlinear Riemann problem 
is solved at each interface between grid zones, and the wave so- 
lutions are used directly to update the mean dynamic variables 
in each zone. In regular zones not including a detected sub- 
shock the gas pressure is updated from the modified entropy 
as P g = S /o 7s_1 , where S is updated by equation (5) outside the 
gas shock (subshock). Updating the gas energy from equation 
(5) instead of equation (3) minimizes errors introduce to the 
gas pressure due to dominant kinetic energy or P c gradients in 
the CR shock precursor. For regular zones including a detected 
subshock, gas pressure is computed in the more conventional 
manner from the updated internal energy found by equation (3). 
The latter is found by subtracting the kinetic energy estimated 
from equation (1) and (2) from the updated total energy (omit- 
ting the d(P c u)/dx term across the subshock jump, since by 
construction the diffusive CRs interact with the gas only over 
distances exceeding the physical shock thickness). 

The employed shock tracking method locates the exact posi- 
tion of the shock within a regular grid zone using the shock 
speed obtained from the nonlinear Riemann solver and then 
adds a new zone boundary there. This subdivides a uniform 
zone of the underlying regular grid into two sub-zones. In the 
next time step, this zone boundary (shock front) is moved to a 
new location using the Riemann solutions at the current shock 
location and the waves are propagated onto the new set of grid 
zones (LeVeque and Shyue 1995). Since the new grid locates 
the shock precisely at an irregular zone boundary, the shock re- 



mains as an exact discontinuity without smearing. The jump 
conditions across the discontinuous subshock are realized ex- 
actly, because the wave solutions from the nonlinear Riemann 
solver (i.e., speeds of waves and jumps associated with three 
wave modes) are used directly to update the dynamic variables 
in each zone. Physically, the shock is not a discontinuity, of 
course, but a relatively thin layer in which entropy is generated. 
Since CR shocks are collisionless, this region is difficult to de- 
fine precisely. Nonetheless, a basic concept behind diffusive 
shock acceleration theory is that we can identify this region and 
that it is thinner than the interaction scale for even the lowest 
energy CRs, given roughly by L>diffO?inj)- The need to satisfy 
this constraint is the main driver for incorporating the shock 
tracking scheme into the CRASH code, since it minimizes the 
number of grid zones containing the subshock. 

In the CR-dynamic step, dP c / dx terms are added to the gas 
momentum conservation equation (2). The gas total energy, 
e g is updated by adding the new, CR-corrected, kinetic energy 
to the thermal energy. At the subshock (only) the gas energy 
must also be reduced according to the term L in equation (11) 
to account for energy removed from the thermal population as 
particles are converted into CRs (i.e., as they are injected). The 
quantity L represents the energy per unit volume per unit time 
transferred from gas to CRs in the immediate postshock region. 
Physically this should occur over a distance corresponding to 
the scattering lengths of injected particles, or roughly Ddiff(/>mj)- 
In practice we spread this exchange over three postshock zones 
of the finest refined grid (see the following paragraph), which is 
slightly broader, but reduces numerical difficulties encountered 
if the energy is extracted too abruptly. This step is analogous 
to the inclusion of radiative energy losses in shocks, except that 
here the energy extraction is limited to very close postshock 
distances. 

In order to follow accurately the evolution of a CR modified 
shock, it is necessary to resolve the precursor structure imme- 
diately upstream of the subshock and, at the same time, to solve 
correctly the diffusion of the low energy particles near the in- 
jection pool. Both of those have scales slightly greater than 
Ddis(Pmj)- While the full extent of the precursor increases with 
the diffusion length of the maximum accelerated CR momen- 
tum, the dominant scale length of the precursor as determined 
by the gas pressure distribution is similar to an averaged diffu- 
sion length of the particle populations with the greatest contri- 
bution to the CR pressure. As a result of these complications, 
a large dynamic range of resolved scales from the shock thick- 
ness to the precursor scale is required for CR shock simulations. 
To solve this problem generally, in order to allow eventually 
for arbitrary and self-consistent diffusion models, we have suc- 
cessfully combined a powerful "Adaptive Mesh Refinement" 
(AMR) technique (Berger and LeVeque 1998) and aforemen- 
tioned "shock tracking" technique and implemented them into 
a hydro/CR code based on the wave-propagation method (Kang 
et al. 2001). The AMR technique allows us to "zoom in" inside 
the precursor structure using a hierarchy of small, refined grid 
levels applied around the subshock. The shock tracking tech- 
nique follows hydrodynamical shocks within regular zones and 
maintains them as true discontinuities, thus allowing us to re- 
fine the region around the gas subshock at an arbitrarily fine 
level. The result is an enormous savings in both computational 
time and data storage over what would be required to solve the 
problem using more traditional methods on a single fine grid. 

Fig. 2 shows the evolution of the gas density within the re- 
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finement region on the l g = 1,3,5,7 level grids for a Mach 100 
shock (see §4.4 for detailed model parameters). There are 200 
zones in each refined grid level and 2 x 10 4 zones in the base 
grid. The size of the refinement region is reduced by a factor of 
two between adjacent grid levels, so the spatial resolution is re- 
fined by the same factor. Fig. 2 demonstrates that the refinement 
region shrinks by a factor of 2 2 for A/ g = 2 and moves along 
with the shock. We also plotted the values at 200 zones with 
filled circles at t = 20 in order to demonstrate that the subshock 
is resolved as a discontinuous jump (see §4.3 for the definition 
of time variable f). 

4.2. Bohm Diffusion Model 

The spatial diffusion coefficient can be expressed in terms 
of a mean scattering length, A, as K(x,p) = h\v, where v is 
the particle velocity. The scattering length, A, and thus n(x,p), 
should be in principle determined by the intensity of resonantly 
interacting Alfven waves. For example, the Bohm diffusion 
model represents a saturated wave spectrum and gives the min- 
imum diffusion coefficient as k b = l/3r g w, when the parti- 
cles scatter within one gyration radius (r g ) due to completely 
random scatterings off the self-generated waves. This gives 
Kb °c p 2 /(p 2 + 1) 1,/2 . We consider here only the proton CR 
component. As stated before, we express particle momenta in 
units of m pc . In order to model amplification of self-generated 
turbulent waves due to compression of the perpendicular com- 
ponent of the magnetic field, the spatial dependence of the dif- 
fusion is modeled as 

K(x,p) = K(p)(-^-), (14) 

pyx) 

where po is the upstream gas density. This form is also required 
to prevent the acoustic instability of the precursor (Drury and 
Falle 1986; Kang, Jones and Ryu 1992). 

Here we define the acceleration time scale for a particle to 
reach momentum p (e.g., Drury 1983) as 

Taccip)= , P ,, = — ^— ( — + — ) ~ \kb(p) , (15) 

<dp/dt> u\ — U2 mi M2 Mj 

where the subscripts, 1 and 2, designate the upstream and down- 
stream conditions, respectively. The last expression is approx- 
imated for the velocity jump for a strong gas shock. For the 
Bohm diffusion coefficient the acceleration time of relativis- 
tic particles increases linearly with the momentum, so both 
the required simulation time and the simulation box size in- 
crease linearly with the highest momentum to be accelerated. 
Hence it becomes extremely costly to extend the CR spectrum 
to p m ia. >• 1 even with our very efficient CRASH code. 

4.3. Model Parameters 

A CR modified shock has more complex structure than a dis- 
continuous jump. A precursor develops and the subshock weak- 
ens as the flow is modified by the CR pressure. We will follow 
Berezhko and Ellison (1999) to denote the values at the far up- 
stream as uq and po, the values immediately upstream of the 
subshock as u\ and pi, and the values immediately downstream 
of the subshock as 112 and p2- We use u s as the subshock ve- 
locity relative to the immediate pre-subshock flow, and Ud as 
the immediate downstream flow velocity in the rest frame of 
the subshock. The gas flow is compressed adiabatically in the 
precursor, so the gas pressure immediately upstream of the sub- 
shock is given by P g;1 = P g fl(pi //>o) 7s - 



The evolution of the CR modified shock depends on four pa- 
rameters: the gas adiabatic index, j g ; the initial Mach num- 
ber of the gas shock, Mo = «o|/c s ,o; P = u /c; and the dif- 
fusion coefficient, n(p). Here uq is the upstream flow veloc- 
ity relative to the initial, unmodified shock, and u Q is the ve- 
locity normalization constant, which is set to be \uq\. We as- 
sume 7 g = 5/3. For all simulations we present here u Q =3000 
kms" 1 , so /3= 10" 2 . We assume a Bohm-like diffusion coeffi- 
cient, k b = n p 2 /(p 2 + l) 1 / 2 . The length and the time variables 
scale with k through the diffusion length and time defined as 
'"o = k / u D and t = k / m 2 , respectively. Thus we do not need to 
choose a specific value for long as we adopt r a and t a as 

normalization constants for length and time scales, respectively. 

We consider a wide range of initial Mach numbers, Mq = 2 - 
200 for the initial shock jump by adjusting the initial preshock 
gas pressure by 

P g:0 = 0.75/(1.25M 2 -0.25). (16) 
The far upstream density and flow velocity are, in code units, 
fixed at po = 1, wo = — 1. Thus the preshock gas is colder for 
stronger shocks, while the gas density and the shock velocity 
are the same for all models. Then the initial jump conditions for 
downstream region in the rest frame of the shock are given as: 
92 = r, U2=—\/r, P & 2 = 0.75 in downstream region, where the 
compression ratio is r = 4M 2 /(Mq + 3). Since the subshock ve- 
locity changes as the subshock weakens in response to growing 
upstream CR pressure, the subshock drifts in the initial shock 
rest frame. In order to keep the subshock close to the center of 
the simulation box a small velocity shift (M s hif t ), which is exper- 
imentally determined, is added occasionally to the rest frame of 
the initial shock. There are no pre-existing CRs in the simulated 
flows, so the CRs are introduced to the system at the shock via 
thermal leakage as the system evolves. 

Although the theoretically preferred values of e lie between 
0.3 and 0.4 for strong shocks (Malkov 1998), such values lead 
to very efficient initial injection for shocks of Mq ^ 100. Effi- 
cient injection at the start of the simulation results in the rapid 
surges of P c and the precursor at the early stage, which, in fact, 
are too fast for the code to follow accurately. During this early 
development the gradients in the precursor are extremely steep. 
Consequently, the flow structure on the base grid, which is cal- 
culated from the average of 2' m!,x zones of the finest grid, devel- 
ops large discrete jumps. The Riemann solver for such flows 
can fail to converge. For e = 0.25-0.3 we can use up to l mux = 5 
levels of refined grids for M < 100 shocks without failures of 
the Riemann solver. For e = 0.17-0.2, l max = 7 can be used for 
Mq < 200 shocks. So we consider two sets of models depend- 
ing on the injection rate, that is, models with e = 0.25-0.3 and 
models with e = 0.17-0.2. While these parameter values were 
applied for practical reasons, the thermal leakage is expected 
to be more difficult in oblique shocks, so these runs can offer 
some insights into the importance of that effect. 

For the injection models with e = 0.25-0.3, the simulations 
were carried out on a base grid with Axo = 3.2 x 10" 3 using 
'max = 5 additional grid levels, so Axs = 10" 4 on the finest grid. 
Here spatial resolutions are given in units of r a . The simu- 
lated space is x = [-12.8, 12.8] and N = 8000 zones are used 
on the base grid. We have repeated one of the models (M = 30 
case) with l m - dx = 4 and found that the numerical results are well 
converged within 1 % in terms of CR pressure at this spatial 
resolution. The number of refined zones around the shock is 
N r f = 200 on the base grid and so there are 2N r f = 400 zones on 
each refined level. The length of the refined region at the base 
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grid is 0.64, so 1 /40 of the entire simulated space on the base 
grid is refined. In these cases we integrate the simulation until 
t = t/t = 20, so that the maximum momentum achieved by the 
end of simulation is of order of p mnx ~ 4, above which the CR 
distribution function decreases exponentially. 

For the injection models with e = 0.17-0.2, the base grid 
spacing is Ax = 1.28 x 10" 2 with Z max = 7, so Ax 7 = 10" 4 on 
the finest grid. The simulated space is x= [-128., 128.] and 
N = 20000 zones are used at the base grid. For these models 
N r f = 100 on the base grid and 2N r f = 200 at each refined level. 
The length of the refined region at the base grid is 1.28, so only 
1 /200 of the entire simulated space on the base grid is refined. 
We integrate these injection models until t = 100-300. For all 
models we use 230 logarithmic momentum zones in the interval 
\og(p/m p c)=[\ogpQ 1 \ogp x \ = [-3.0,+3.0]. 

Gasdynamic variables are assumed to be continuous at both 
left and right boundaries of the simulation box. For the CR 
distribution function a continuous boundary is assumed for the 
advection step and a no-flux boundary condition is adopted for 
the spatial diffusion step. Either below po or above p\ ; g(p) = 
is assumed. For clarity we distinguish here between the highest 
momentum included in the simulations, p\ and the value of the 
momentum /? max characterizing the effective upper cutoff in the 
actual distribution, g(p). In practical terms there generally is 
a clear momentum above which g(p) drops sharply, so that we 
identify as Pmax- The maximum momentum p m - dx will gener- 
ally evolve towards, p\, but in all the simulations reported here 
remains well below it at the terminal time. 

4.4. Injection and Acceleration Efficiencies 

The fraction of particles that have been swept through the 
shock after the time t, and then injected into the CR distribu- 
tion is estimated by 

/ dx / 47t/cr0? , x, t )/? 2 dp 



(17) 



/ n «o( f ')dt' 

where /cr is the CR distribution function, «o is the particle 
number density, and u' Q (t) is the instantaneous shock velocity 
relative to the far upstream flow. As the CR pressure becomes 
dominant for strongly CR mediated shocks, the effective adia- 
batic index of the flow becomes smaller than 5/3. This com- 
bines with other factors discussed below that allow the total 
shock compression to increase and causes the subshock to slow 
down with respect to the upstream gas. So the instantaneous 
shock velocity decreases as the CR pressure increases. One 
should also note that £(/) represents a "moving time-averaged" 
injection efficiency rather than a current instantaneous value. 

As a measure of acceleration efficiency, we define the "CR 
energy ratio"; namely the ratio of the total CR energy within 
the simulation box to the kinetic energy in the initial shock rest 
frame that has entered the simulation box from far upstream, 

/ dxE CR (x,t) 



$(*) = ■ 



0.5po|«o| 3 f 



(18) 



The total energy in the simulation box is conserved in the initial 
shock rest frame, since the total energy flux entering the sim- 
ulation from far upstream, 0.5/£>oKo + 7g^g,oWo/(7g- 1), is bal- 
anced by the same flux exiting the grid far downstream. We 
intentionally placed the boundaries well away from the shock 
to avoid boundary issues. Since our shock models have the 
same upstream density and velocity, but different gas pressure 
depending on Mo, we use the kinetic energy flux rather than the 



total energy flux to normalize the "CR energy ratio". For strong 
shocks the total energy flux becomes the same as the kinetic en- 
ergy flux. 

Since the shock slows down due to nonlinear modification, 
the kinetic energy flux in the instantaneous shock rest frame 
also decreases. So we also calculate the ratio of the total CR en- 
ergy in the simulation box to the kinetic energy defined in the 
instantaneous shock frame which has entered the shock from 
far upstream, 

JdxE CR (x,t) 
* W /0.5p « (0 3 dt- 

5. SIMULATION RESULTS 

5.1. Models with e = 0.25 -0.3 
5.1.1. Shock Structure Modification 

Figs. 3 and 4 show time evolution of the CR modified shock 
structure for Mo = 5 and Mo = 30 shocks, respectively, with 
e = 0.3. As CRs are injected and accelerated, the CR pressure 
increases, the precursor grows, and the subshock slows down 
and weakens. Nonlinear modification to the shock structure 
is much stronger and the subshock drift is greater for higher 
Mach number models. At the simulation ends the subshock 
drifts to the left with warm ~ -0.14 for the Mo = 30 shock and 
with Mdrift ~ -0.038 for the Mo = 5 shock. Since the initial sim- 
ulation frame was set up with the w S hitt = 0.01 relative to the ini- 
tial shock rest frame, the subshock drifts in the initial shock rest 
frame with u = Mdrift _ 0.01. Thus for strong shocks the subshock 
velocity relative to far upstream flow becomes u' ~ 0.85 instead 
of | mo | = 1, approximately independent of Mo. For low Mach 
number shocks, however, this velocity difference is smaller and 
increases with Mo (see also Fig. 12 below). 

Due to strong injection, the CR pressure increases and the 
modification to the flow structure proceeds very quickly, be- 
fore f < 0.1. After the initial quick adjustment, the CR pressure 
at the shock reaches approximate time-asymptotic values when 
the fresh injection and acceleration are balanced with advection 
and spreading of high energy particles due to strong diffusion. 
This balance seems to occur earlier in time for larger values 
of e for models of a given Mach number (see also Fig. 5 be- 
low). For a given value of e, it occurs gradually earlier with in- 
creasing initial Mach number, but at a similar time in the strong 
shock limit. We will show below that the time-asymptotic val- 
ues of P c are related to the subshock evolution, as well. After 
P c 2 becomes quasi-steady, the overall shock structure evolves 
approximately in a self-similar way, as shown in Figs. 3-4. 
This behavior can be understood from the fact that P C 2 at the 
shock stays constant and the characteristic scale length of P c 
and precursor scale length increase linearly with time, since the 
diffusion length of the highest momentum, £>diff(/?max) oc t for 
Pmax ^ 1 and the advection length also scales with time. 

The evolution of the density distribution in Fig. 4 includes 
the formation of the postshock "density spike" noted in sev- 
eral previous time dependent studies (e.g., Jones and Kang 
1990). The density spike forms because the weakening of the 
subshock lags behind formation of the precursor during initial 
shock modification. For strong shocks this can produce a sub- 
stantial, temporary "overshoot" in the total compression (see 
also Fig. 9). That feature eventually detaches from the steady 
shock and is carried downstream from it. 

Evolution of subshock and precursor properties can be seen 
in Fig. 5, where the subshock velocity, u s and Mach number, 
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M s , as well as the pre-subshock gas density, p\ j po, and the im- 
mediate postshock gas density, pi/po, are plotted against time 
for the models of Mo = 2- 30 and e = 0.3. As mentioned be- 
fore, the modification to the flow structure occurs mostly before 
t < 0.1, and so the subshock velocity decreases and the precur- 
sor compression increases at the same time. After t > 0.1, the 
subshock Mach number remains nearly constant at M s <~ 2 - 3 
for these shock models, and the subshock persists. Thus, a 
completely smooth transition never develops. The ratio p 2 / Po 
is greater for higher Mach number shocks and in fact follows 
roughly M{J 6 scaling relation (see §5.3 for further discussion). 
In our simulations, the compression continues to grow slowly 
with time, especially for strong shocks, even after the overall 
shock structure measured from far upstream through the sub- 
shock appears mostly to have reached a "quasi-steady state". 

Obviously the compression factor is much higher than what 
is expected in a steady-state, energy conserving shock jump 
for relativistic gas (i.e., p 2 / po = 7 with j c = 4/3). Similar re- 
sults to ours have been found previously in spherically expand- 
ing CR dominated shocks (Berezhko et al. 1995) and in steady 
state planar shocks with strongly momentum dependent diffu- 
sion when provision is made for energy losses carried by es- 
caping high momentum particles (Malkov 1997; Berezhko and 
Ellison 1999; Malkov and Drury 2001). While our shocks are 
not yet truly steady, they have all evolved well past the tran- 
sient stage associated with the density spike mentioned above. 
It is very unlikely that the density jumps across our simulated 
shocks would decrease were they allowed to evolve over much 
longer times. As noted, compression appears still to be increas- 
ing. On the other hand, the spatial and momentum boundaries 
in our simulations are such that energy losses for the entire sys- 
tem can be ignored, so they might be considered to be conser- 
vative. 

The key feature in common to the CR shock solutions with 
non-transient large density jumps is that they all involve sig- 
nificant energy transport away from the transition itself. Drury, 
Volk and Berezhko ( 1 995) pointed out that subshock smoothing 
in spherical systems was inhibited by the fact that CRs diffuse 
into an ever-increasing volume, thus extracting energy from the 
shock itself. Steady planar shocks with either a spatial or a 
momentum "escape" boundary that allows significant energy to 
leave the shock structure encounter the same effect. Once the 
momentum distributions in our CR populations near the shock 
evolved so that p max ~ p\ they would be influenced by similar 
"escape" losses. In the meantime they appear to be evolving 
towards properties that are consistent with those losses. At the 
"intermediate ages" we see in our shock evolutionary histories, 
energy is not leaving the system, but is being transported out of 
the region of the subshock by diffusion of the highest energy 
CRs. The quasi-steady structure is growing in a self-similar 
manner, in fact, reflecting the linear increase in D^(p m . dx ) with 
time. Since advection extends the downstream region in a sim- 
ilar way, the immediate postshock properties can remain essen- 
tially constant even while the global structure continues to grow 
in phase space. 

Regardless of the quasi-steadiness of the postshock CR pres- 
sure and the self- similarity of the flow structure, the CR dis- 
tribution function continues to become harder and extends to 
ever higher momenta until scattering becomes too weak to 
contain the particles. Thus, the adiabatic index of the CRs, 
7 C = 1 +P C /E C , decreases with time at the shock and tends to- 
wards 4/3 as the CR pressure and energy become dominated 



by relativistic particles. This evolution in "f c is more dramatic 
away from the subshock and especially upstream, since spa- 
tial diffusion is biased towards more energetic particles. This 
secular evolution of j c affects the precursor compression and 
causes the subshock properties to evolve slowly, especially for 
strong shocks, as can be seen in Fig. 5. Although we adopted 
the modified entropy equation to follow the precursor adiabatic 
compression accurately, the subshock properties are not numer- 
ically steady but somewhat noisy. This unsteady behavior of the 
subshock seems to be due to the CR pressure gradient, dP c /dx, 
which is separately added to the momentum conservation equa- 
tion. Within the precursor this pressure force is large and com- 
presses the gas adiabatically, so it determines the position and 
the properties of the subshock. Thus, inside the precursor the 
accuracy in calculation of P c distribution would control the sub- 
shock calculation rather than the accuracy in calculation of the 
gas adiabatic compression. If one used a conventional code 
with the total energy equation alone, these quantities would be 
much noisier than what is seen in Fig. 5. 

5.1.2. Injection and Acceleration Efficiencies 

Fig. 6 shows how the CR energy ratio, $, the CR pressure 
at the shock normalized to the far upstream ram pressure in 
the instantaneous shock frame, Pc.i/ Pou' Q (t) 2 , and the "time- 
averaged" injection efficiencies, £, evolve for shocks with dif- 
ferent Mach numbers when e = 0.3 (left three panels) or e = 0.25 
(right three panels) is adopted. As mentioned in the previous 
section, for all Mach numbers the postshock P c increases until a 
balance between injection/acceleration and advection/diffusion 
of CRs is achieved, and then stays at a steady value after- 
wards. The time-asymptotic value of the CR pressure becomes, 
Pel/ Pou'oit) 2 ~ 0.8, forMo = 30 with e = 0.3. As seen in the pre- 
vious section, the subshock slows down, and the subshock ve- 
locity relative to the far upstream flow becomes u' Q ~ 0.85 |«o| 
for Mq = 30. Since the subshock has reached a "quasi-steady 
state" in our simulations, we can apply the mass and momen- 
tum flux conservation in the subshock rest frame as 

p u' = p\u s = p 2 u d , (20) 

P2U 2 d + Pg,2 + Pc,2 = P0UQ+Pgfl- (21) 

The postshock gas pressure is given by the jump condition, 
P gi 2 = (1.25MJ -0.25)P g .i, so it can be approximated as P g2 ~ 
0.15p\u 2 for 2 < M s < 4 for models considered here (see Fig. 
5). Hence for large M the postshock CR pressure can be 
approximated as P c , 2 /pou' ( 2 ~ 1 -0.75(m j /m )-(mj/mo). Since 

u'q « 0.85|«o| and u s sa 0.17|ko| forMo = 30 shock, P C}2 / Pou'q ~ 
0.8, which is consistent with what we found in our simulations. 

The CR energy ratio, $, increases with time, but asymp- 
totes to a constant value, once P Ct2 has reached a quasi-steady 
value. As mentioned in the previous section, this results from 
the "self-similar" behavior of the P c distribution. Since both 
the numerator and denominator in equation (18) increase ap- 
proximately linearly with time in the case of the self-similar 
evolution of E c , $ becomes steady when P c2 becomes con- 
stant. Time-asymptotic values of $ increase with the initial 
Mach number M and $ w 0.53 for M = 30 at the terminal 
time (see also Fig. 12). On the other hand, time asymptotic 
values of $' become $' w $ • (m /|m |) 3 ~ 0.86 for the M Q = 30 
shock model. 

The initial injection rate before the weakening of the sub- 
shock is roughly independent of Mo for strong shocks, as we 
expect from the discussion in §3.3. Initially it briefly increases 
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to £ ~ 0.01 for strong shocks until t ~ T acc (/?; n j), the accelera- 
tion time scale of the injection momenta. The initial brief rise in 
£ is a numerical artifact occurring in the interval when the ini- 
tially pure Maxwellian distribution adjusts to form a nonther- 
mal tail in the momentum range where injection takes place. 
Afterwards it decreases as the subshock weakens, since the 
injection is less efficient at low Mach number shocks. After 
the initial adjustment period the injection rate decreases as a 
power-law form with time, £ ~ £ f~ a for 0.05 < i < 20, where 
0.15 a 0.4 with larger values for higher Mq. Although 
the subshock Mach number decreases to similar values after 
/ = 0.1 for all models (see Fig. 5), the injection rate decreases 
faster for higher Mq models. This can be explained as follows. 
The injection momentum p m ^ after t = 0.1 is smaller for higher 
Mq, since the subshock velocity is reduced further in response 
to stronger structure modification. Diffusive particle flux up- 
stream, which is realized numerically in our code with a given 
subshock thickness (Axs), is reduced if AuffCPinj) < AX5. The 
resulting injection rate is smaller for higher Mo models in which 
the injection momentum is lower. Thus the strong decrease of 
£ for high Mq models below the level of £ for Mq = 3 model 
should be regarded as a numerical artifact due to a fixed shock 
thickness. 

In order to explore the dependence of our injection model 
on the parameter e, we also considered the same shock models 
with e = 0.25 and show $, P c ,2 and £ in the right three panels of 
Fig. 6. Because of less efficient injection compared to e = 0.3 
cases, the initial increases of P c 2 and $ are delayed somewhat 
and their time-asymptotic values are reduced a little, as well. 
Before t <~ 1, the CR pressure is dominated by non-relativistic 
particles and the details of the injection history are important to 
the shock structures. Afterwards, the relativistic particles dom- 
inate and the results become less sensitive to the values of e and 
the details of the injection process. For example, the injection 
rate £ is 2-3 times lower with e = 0.25 than with e = 0.3 for the 
Mq = 30 shock, but the time-asymptotic value of $ is smaller 
only by 2%. In summary, the CR energy ratio, $, depends sen- 
sitively on Mq for weak shocks, then becomes independent on 
Mq for strong shocks, but depends only weakly on the injection 
rate (e). 

5.1.3. The Particle Distribution 

Time evolution of the distribution function g(p) = f(p)p A at 
the shock is shown in Fig. 7 (Mq = 5) and Fig. 8 (Mq = 30) 
along with the changes in the transparency function r esc . The 
Maxwellian distribution shifts to lower momenta primarily due 
to weakening of the gas subshock in response to flow deceler- 
ation in the precursor, but also due to the loss of thermal en- 
ergy to the injection process. As the subshock weakens the 
postshock flow speed in the subshock rest frame, mj, decreases, 
so the transparency function also shifts to lower momenta. As 
shown in Fig. 1, the ratio of w; n j/wth increases as the subshock 
weakens, leading to decreases in the injection rate. For mo- 
menta just above the injection pool, the distribution function 
changes smoothly from a Maxwellian distribution to an approx- 
imate power-law whose index is close to the test-particle slope 
for the subshock. However, we note that the particle distribu- 
tion at lower energy does not evolve as a simple power-law, 
but, instead, it depends on the time-dependent injection his- 
tory, since the subshock properties and r esc (Mj) vary with time 
and their numerical solutions in our simulations contain low- 
amplitude noise. Those "historical details" become smeared 



out at higher momenta. While a fraction of particles injected 
earlier continues to be accelerated to higher momenta, so that 
Pmnx(t) increases, the amplitude of g(p) at the shock and at a 
given momentum decreases with time for p < p mix (t) due to 
diffusion. 

These results can be compared with those of Gieseler et al. 
(2000) in which a hydrodynamic code based on the Total Vari- 
ation Diminishing (TVD) scheme was used. In Gieseler et 
al. (2000) the diffusion coefficient was 100 times larger than 
the value used in the present study, so in a comparison the 
length and time scales should be scaled by that amount. The 
main difference between the current CR/AMR code and the 
CR/TVD code is that the subshock is maintained as a discontin- 
uous jump in the CR/AMR code. Numerical gradients defined 
at the subshock, therefore, span a single regular zone. Given 
the same nominal spatial resolution, the shock is more sharply 
resolved in the CR/AMR simulations. This affects the diffusion 
of suprathermal particles near the shock and so leads to higher 
injection rate. 

5.2. Models with e = 0.17-0.2 

In order to extend the calculation to much longer simulation 
times, so that much more energetic particles will be accelerated, 
we need to include a larger computational domain because of 
the fast spatial diffusion of the highest energy particles. Thus, 
the computational time for such calculations becomes very de- 
manding. We increase the refinement level to / max = 7 and con- 
sider a computation domain 10 times larger than that of the 
stronger injection calculations discussed in §5.1. This forces 
us to choose smaller values of e = 0.17-0.2 (weaker injection 
rate) to avoid the initialization numerical problem mentioned in 
§4.3. For Mq = 5 - 100 e = 0.2 is used, while for M =150 and 
200 e = 0.17 is used for t < 1 then increased to e = 0.2 for t > 1. 

5.2.1. Shock Structure Modification 

Fig. 9 shows the evolution of the shock structure for a 
Mq = 200 shock at much later times up to t = 140. Overall shock 
structures reach approximate time-asymptotic states after the 
initial quick adjustment, and then follow "self-similar" evolu- 
tions afterwards as in the the Mq = 30 shock with e = 0.3 shown 
in Fig. 4. The subshock slows down and drifts with the velocity 
Mdrift ~ -0.025 after the shock structure evolution became quasi- 
steady for the strong shock models and the starting velocity 
shift relative to the initial shock rest frame was M s hif t = 0. 145. So 
the time-asymptotic velocity of the subshock relative to far up- 
stream becomes u' () s=s 0.83 for e = 0.2 for strong shocks (see also 
Fig. 12). This shock speed was u' Q m 0.85 for Mq = 30 shock 
model with e = 0.3. The subshock Mach number decreases with 
time, but reaches quasi-steady values after t = 100. The values 
of M s at t = 100 seem to scale with the initial shock Mach num- 
ber Mq roughly as M s - 2.9MJJ 13 . Although the postshock CR 
pressure and the subshock speed become roughly quasi-steady 
after the initial fast modification stage, the precursor compres- 
sion continues to increase up to the terminal time and the rate 
of this increase is higher for higher Mq models. We see that 
7e is the largest near the shock where the injected suprathermal 
particles are dominant. But it becomes smaller away from the 
shock, because relativistic particles dominates there. The spa- 
tial distribution of P c peaks sharply from upstream at the sub- 
shock on account of the small diffusion length scale reflecting 
the dominance of low energy particles near the subshock. 
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5.2.2. Injection and Acceleration Efficiencies 

The evolution of $,P C , 2 , and £ for M = 5-200 with e = 0.2 is 
shown in the left three panels of Fig. 10. Note again that we set 
e = 0. 17 when t < 1 and e = 0.2 when t > 1 for M = 150 and 200 
in order to avoid the initial convergence failure of the Riemann 
solver. So the injection is slower in these two models, result- 
ing in a retarded evolution of all related quantities until t > 1 . In 
the evolved shocks the CR pressure at the subshock is higher for 
larger values of Mq and becomes P c 2 « 0.9pqu' q 2 f° r Mq <; 50. 
For high initial Mach-number shocks (Mq > 50) the CR energy 
ratio approaches $ « 0.6, resulting in <&' w $ • (|mo|/« ) 3 ~ 1> 
since the approximate time-asymptotic value of m /|«o| is 0.83 
for strong shocks. This high efficiency is possible, because CRs 
have absorbed some of the thermal and kinetic energy from the 
initial shock structure inside the simulation box, in addition to 
the kinetic energy that has entered the shock during the simula- 
tion, / 0.5pQU 3 dt. We note that the "undulating" features in the 
time evolution of P C 2 seem to be numerical artifacts, and not 
real. Unlike other spatially averaged or integrated quantities, 
the plotted P c ,2 is sampled exactly at the subshock (i.e., from 
one zone) whose properties show small, noisy variations in 
time. These features seem in particular to correspond to times 
when the subshock crosses a regular zone boundary and are 
most prominent for the injection models with e = 0.17 - 0.2. 

Once again the initial injection rate before the weakening of 
the subshock is independent of Mq for strong shocks, and it is 
about 4 times smaller than that of e = 0.25 models and about 
8 times smaller than that of e = 0.3 models. For 1 < t < 200, 
the time averaged injection rate can be fitted as a power-law, 
£ <~ 2.5 x 10" 3 f -0 - 4 , for strong Mach shocks. As discussed in 
§5.1.2, the injection rate decreases with time mainly because 
the subshock weakens due to CR modification, but partly be- 
cause the diffusive flux of injected particles crossing the shock 
decreases for smaller injection momentum. If the effective nu- 
merical subshock thickness (spatial resolution) were set up to 
decrease as the postshock temperature decreased, we expect 
the injection rate for high Mq shocks would level off to the 
level similar to that for low Mq shocks, that is, £ ~ 10" 3 4 found 
for the Mq = 5 case. Below we will list this for our estimated 
asymptotic rate. Considering the lack of our understanding of 
the detail physics of the shock transition, this should be consid- 
ered a rough estimate. 

The same set of quantities are shown in the right three panels 
of Fig. 10 for Mq = 30 for e = 0.2, 0.23, and 0.25, in order to 
study the effects of different injection rates. Obviously larger 
values of e (i.e., more leakage) lead to higher injection rate, £, 
and so acceleration of CRs proceeds faster. For all values of 
e considered, however, time-asymptotic values of P c ^/ ' pq(uq/) 2 
~ 0.8 and $ ~ 0.6, even though the injection histories are dif- 
ferent. This implies that the long term behavior of $ is mostly 
dependent on the initial shock Mach number, but only weakly 
dependent on the details of the injection rate, so long as it is 
in the range seen here. If the injection rate falls below some 
critical range, other studies have emphasized that shock modi- 
fication can be minimal, so that "test particle" solutions apply 
(e.g., Berezhko et al. 1995; Malkov, Diamond and Volk 2000). 
It seems mostly likely that all the shocks simulated here will 
remain strongly modified, so correspond to "high injection effi- 
ciency", although we cannot yet judge if they would eventually 
approach some kind of self-organized critical condition such as 
that suggested by Malkov, Diamond and Volk (2000). Evaluat- 
ing that possible eventuality will likely require self-consistent 



treatments of the coupling between evolution of CRs and the 
diffusion coefficient. That is beyond the scope of the present 
work. 

5.2.3. The Particle Distribution 

The evolution of the CR distribution function for a Mq = 100 
shock is shown in Fig. 1 1 . The general behavior is the same as 
the early evolution for the Mq = 30 shock with e = 0.3 shown 
in Fig. 8. In particular both the thermal distribution and the 
transparency function shift to lower momenta, the function g(p) 
shows the characteristic "concave upwards" curves reflecting 
modified shock structure. The CR spectrum extends to higher 
momenta with time, but the amplitude at intermediate momenta 
decreases with time. The slope of the CR spectrum, defined as 
q = -(dlnf /dlnp), ranges over 4.1 < q < 4.4 near p m j, reflect- 
ing r su b > 3 then decreases with the particle momentum and 
converges for strong shocks to q ~ 3. 1 just below /? max - 

5.3. Time Asymptotic Behaviors 

Although our CRASH code is much more efficient than con- 
ventional codes used for this problem, we have so far integrated 
our models only until the maximum momentum reaches about 
~ 40m p c (t = 200), since computational requirements of using a 
Bohm type diffusion model are still substantial. However, some 
characteristics of the shock seem to have reached quasi-steady 
states well before the end times of our simulations. So we can 
make an attempt to project the long-term, time-asymptotic be- 
haviors of these CR modified shocks from our results. For 
strong shocks of Mq > 30 the CR injection via thermal leak- 
age is very efficient, and so the CR pressure quickly domi- 
nates over the gas pressure. The postshock CR pressure reaches 
approximate time-asymptotic values by the time the injection 
and acceleration become balanced by the advection and diffu- 
sion of CRs away from the shocks. As compression increases 
through the precursor, the subshock slows down, and the sub- 
shock velocity relative to the far upstream flow approaches 
u' ~ (0.83 — 0.85)|«o| for strong shocks with the inverse wave 
amplitude parameter considered here, e = 0.2-0.3. As a result, 
the postshock CR pressure asymptotes to a quasi-steady value 
of P c .2 ~ 0.9pqu'q 2 ~ 0.63po«o f° r strong shocks. Approximate 
estimates for time-asymptotic values of the CR energy ratio, 
are calculated at t = 20 for models with e = 0.3 and 0.25, and 
at t = 100 for models with e = 0.2. They are plotted as a func- 
tion of the initial shock Mach number Mq in the right top panel 
of Fig. 12. Also time-asymptotic values of u' / 1 uq | = M' (j /Mq, 
are calculated at t = 100 for models with e = 0.2 and shown 
in the right middle panel of Fig. 12. We can estimate time 
asymptotic values of $' as $ • (\uq\/u' q ) 3 using data in these two 
plots. We also show time-asymptotic values of the CR pressure 
at the shock normalized to the ram pressure of the initial shock, 
Pc,i/ PqWq (solid line), and that normalized by the ram pressure 
of the instantaneous shock, P c .2/pqu' ( 2 (dotted line), which are 
calculated at t = 100 for e = 0.2 shocks in the right bottom panel 
of Fig. 12. 

Some other characteristics of the subshock, on the other 
hand, show slow secular changes, even after the postshock P c 
becomes quasi-steady. As the CR distribution continues to 
harden, for example, the effective adiabatic index of the flow 
decreases with time. Also the precursor spreads further out, the 
compression through the precursor increases and the subshock 
continues to weaken. The entire shock structure broadens lin- 
early with time in a "self-similar" way, because of the strongly 
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momentum dependent diffusion matching the advective spread- 
ing behind the shock. As a result, unless we enforce an escap- 
ing upper momentum boundary (i.e., g(p) = for p > p\, where 
pi is the upper momentum boundary), the CR modified shock 
would not reach a strict steady-state. 

According to kinetic equation simulations of spherical blast 
waves (Berezhko et al. 1995; Drury, Volk and Berezhko 1995), 
and both Monte Carlo and kinetic equation calculations of 
steady planar shocks with strongly momentum dependent dif- 
fusion (Malkov 1997; Berezhko and Ellison 1999), a CR mod- 
ified shock characteristically retains a subshock with the com- 
pression ratio 2.5 < r su b < 4 and does not become a completely 
smooth transition. In the present time dependent planar shock 
simulations with Bohm-like diffusion we also found that the 
subshock remains. The values for the subshock Mach numbers 
in our shocks can be described by the scaling M s ~ 2.9MJJ 13 
at the end of the simulations. Since the subshock compres- 
sion ratio is related to the gas subshock Mach number as r su b = 
4/(1 + 3/M 2 ), this leads to 2.5 < r sub < 4. Steady, adiabatic 
compression through the precursor can be expressed simply as 
p l /p = (M^/M S ) ( T8 +1) / 2 , where M' = M ■ (u' /u ) is the time 
asymptotic value of the shock Mach number with respect to the 
far upstream flow. Then the total compression ratio r tot = P2/ pa 
for a steady CR modified shock with a subshock can be ex- 
pressed in terms of the subshock Mach number and the sub- 
shock compression ratio as 

no, = r mh {M' /M s fl 4 (22) 

(Berezhko and Ellison 1999). Here r sub M; 3/4 =4M^ /4 /(M 2 + 3) 

from the shock jump condition. So, r su bM s ^ 4 = 0.96-1.36 for 
the subshock Mach numbers found in our simulations (i.e., 2 < 
M s < 6). Both Berezhko et al. (1995) and Berezhko and El- 
lison (1999) found that r sub M7 3 ^ 4 ~ 1 for total Mach numbers 
(M' ) considered, so that the total compression ratio scales with 
the total Mach number as r tot <~ (Mq) 3 / 4 . For comparison we 
plotted in the left three panels of Fig. 12 the preshock (pi/po). 
postshock (P2/P0) densities, and subshock compression factor 
(r su b) at t = 20 for models with e = 0.3 and 0.25, and at t = 100 
for models with e = 0.2 as a function of the initial shock Mach 
number Mo. We note once again that the precursor compres- 
sion continues, although very slowly, so these quantities have 
not reached time-asymptotic values at the terminal time of our 
simulations. For our simulated shocks with e = 0.2 the ratio 
Pi I Pa shown in Fig. 12 can be approximated by a power-law, 
p 2 /pa ~ l.5M° 6 for M < 80. But this power-law relation flat- 
tens out for M) > 80. Considering that p 2 has not reached 
a quasi-steady value at t = 20 or t = 100, especially for high 
Mach number shocks, and that the ratio M' /Mq decreases with 
Mo, our empirical relation for r tot vs. M' is likely to steepen 
if we run the simulations for much longer times. As discussed 
in §5.1.1, the scaling relation for the compression ratio results 
simply from the fact that the subshock persists in our simulated 
shocks, which broadens in space self-similarly with time due 
to strongly momentum dependent diffusivity, but are not yet 
stationary. According to Malkov (1997), the total compression 
saturates at the level of (vp mm /p m ^) for M > (vp m!a /pmj) 4/3 , 
where v is his injection parameter. For strong shocks in our 
simulations, i/~4x 10~ 2 and p max /pmj ~2x 10 4 , so the shocks 
considered here are not in the saturation limit. 

6. SUMMARY 



In order to study cosmic-ray modified shocks we have de- 
veloped a new numerical code, CRASH (Cosmic-Ray Amr 
SHock), by implementing a thermal leakage injection scheme 
introduced by Gieseler et al. (2000) into a new hydro/CR code 
with Adaptive Mesh Refinement and Shock Tracking scheme 
by Kang et al. (2001). Our CR injection model is based on the 
"thermal leakage" process at quasi-parallel CR shocks. Injec- 
tion is regulated by the convolution of the population in the high 
energy tail of the Maxwellian velocity distribution of the post- 
shock gas and a transparency function, T esc (p, Ud,e), determined 
by the strength of downstream MHD waves, expressed by a pa- 
rameter, e. The injection rate in our model is then controlled 
largely by the subshock Mach number, M s , since that deter- 
mines the ratio of the postshock flow speed to the breadth of the 
thermal distribution. For weak shocks injection becomes more 
difficult as M s decreases, but is independent of M s for strong 
shocks, when the subshock compression asymptotes. With our 
CRASH code the CR injection and acceleration at astrophysical 
shocks can be simulated numerically even for strongly momen- 
tum dependent spatial diffusion coefficients. 

Using these tools the time evolution of CR shocks has been 
followed with a Bohm type diffusion model. We started from 
the initial conditions for pure gasdynamic shocks of various ini- 
tial Mach numbers (Mo) without any pre-existing CRs. In such 
simulations the CR injection rate is not treated as a fixed free 
parameter, as it has been done traditionally. Instead, once an 
"intelligent estimate" is made of the strength of the trapping 
wave field of the downstream MHD turbulence, the injection is 
followed naturally and self-consistently during nonlinear evo- 
lution of the flows. For strong shocks with initial strengths, 
Mo > 30, a substantial portion of the particles in the tail of the 
Maxwellian distribution have velocities high enough to leak up- 
stream against the wave-particle interactions, so the injection is 
efficient and fast. As the CR pressure increases at the subshock 
and the in-flowing plasma is compressed, however, the sub- 
shock slows down with respect to the far-upstream flow. This 
modification occurs rather promptly and before the postshock 
CR pressure reaches an approximate time-asymptotic value. 
Afterwards the evolution of the shock structure becomes sec- 
ular. For strong shocks, the subshock speed in our simulations 
decreases by about 15-17 %, and the postshock CR pressure ab- 
sorbs up to 60 % of the ram pressure of the initial shock, which 
corresponds to 90 % of the ram pressure of the upstream flow 
in the evolved subshock rest frame at the end of our simula- 
tions. Once the postshock CR pressure becomes constant, the 
shock structure evolves approximately in a "self-similar" way, 
because the scale length of shock broadening increases linearly 
with time. 

The injection rate, defined as the fraction of the particles 
passed through the shock that are accelerated to form the CR 
population, becomes as high as £ ~ 0.01 early in the evolution 
of strong shocks, independent of Mo. As the shock is modi- 
fied, however, the subshock Mach number decreases down to 
M, ~ 2.9Mg 13 , and the injection rate reduces towards the limit- 
ing value corresponding to weak subshocks. Our approximate 
numerical estimate for the limiting value is £ ~ 10~ 3 4 for mod- 
els with e = 0.2. Since the CR spectrum continues to extend 
to higher momenta, and since those highest energy particles 
diffuse rapidly away from the shock when k oc p, energy be- 
gins to be transported away from the shock transition. This 
allows the subshock to persist and the total compression to be- 
come large compared to steady energy-conserving gasdynamic 
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shocks. Also the shock structure evolves towards the high com- 
pression ratios seen in kinetic equation spherical CR shocks and 
steady state planar shocks with energy escape through spatial or 
momentum boundaries. 

Finally, the main conclusions of our time-dependent simula- 
tions of CR modified shocks based on Bohm-like diffusion and 
a physically-based thermal leakage CR injection are: 

1. In the strong shock limit of initial Mach number Mo > 30, 
significant physical processes such as injection and accelera- 
tion become largely independent of the initial shock Mach num- 
ber, and only weakly depend on the postshock wave amplitude 
parameter, e, for values considered here (0.2 < e < 0.3). Ac- 
cording to our thermal leakage model, the overall injection rate 
approaches £ ~ 10~ 3 and the fraction of upstream flow kinetic 
energy in the initial shock rest frame that has been transferred 
to CRs is $ ~ 0.6 for strong shocks. The ratio of the postshock 
CR pressure to the instantaneous ram pressure of the subshock 
with respect to the upstream plasma approaches ~ 90%. On 
the other hand, the injection rate and acceleration efficiency 
are sensitively dependent on Mq for low Mach number shocks 
(M <30). 

2. In our simulations, with no initial CRs around the shock, 
the thermal leakage injection is very efficient initially when the 
subshock is strong, but it becomes much less efficient as the 
subshock weakens due to nonlinear feedback from the CR pres- 
sure. For example, the time-averaged injection rate can be fit- 
ted as a power-law form, £ ~ 2.5 x 10" 3 (f)"°- 4 for 1 < t < 200 
(t = t /to) for strong shocks of Mq J> 30 and the inverse wave- 
amplitude parameter e = 0.2. 

3. Although some particles injected early in the shock 
evolution continue to be accelerated to ever higher energies, 
the immediate postshock CR pressure reaches a quasi-steady 
value when a balance between injection/acceleration and ad- 
vection/diffusion is achieved. The region of quasi-steady post- 
shock properties spreads for Bohm diffusion in an almost self- 
similar fashion because both diffusive and advection rates scale 
linearly with time. This expansion maintains the large precur- 
sor compression at values that can be large compared to what 



one would derive from Rankine-Hugoniot relations for steady 
gas shocks with a relativistic equation of state. 

4. The sum of the compression through the precursor and 
across the subshock calculated near the terminal time, t = 100, 
for models with e = 0.2 can be approximated by P2/P0 ~ 
1.5M ) (5 for M Q < 80, but this p 2 /po-M Q relation flattens for 
Mq > 80. Considering that p 2 continues to increase until the 
terminal time, especially for higher Mach number shocks, and 
that the Mach number of the subshock with respect to the 
far upstream (M' ) decreases for strong shocks with high ini- 
tial Mach numbers (Mo), this is probably consistent with the 

P2/P0 oc M^ A scaling relation found for steady shocks and 
for spherical shocks by several authors (Berezhko et al. 1995; 
Malkov 1997; Berezhko and Ellison 1999). The large density 
compression is possible in the planar shocks without energy es- 
cape from the whole computed system, since the shock struc- 
ture spreads out far upstream and far downstream, due to strong 
momentum dependent diffusivity. 

Although our simulations have been carried out only until 
the highest momentum is p m . dx ~ 40m p c due to severe com- 
putational requirements, these simulation results provide use- 
ful guidance for long-term evolution of CR modified shocks. 
In a future study, we will extend the integration time so that 
Pmax/mpC » 1 can be achieved and the free escape at the up- 
per momentum boundary is applied in order to compare the 
time-dependent simulations with previous studies of steady- 
state shocks (Ellison and Eichler 1985; Malkov 1997; Malkov 
and Drury 2001). 
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APPENDIX 
DISCUSSION OF TEST CALCULATIONS 

The algorithms applied in CRASH are already mostly well-documented and well-tested. The hydrodynamic scheme, for example, 
is a member of the well-known Godunov family and was verified by LeVeque and Shyue (1995). Our Crank-Nicholson scheme 
for solution of the diffusion-convection equation was tested against previous, independent implementations and analytic test particle 
diffusive shock acceleration solutions in Kang and Jones (1991), where we employed the PPM hydrodynamics algorithm. Other tests 
of nonlinear CR shock solutions using several previous hydrodynamical implementations were given for both diffusion-convection 
(e.g., Kang, Jones and Ryu 1992; Kang and Jones 1995) and two-fluid (e.g., Jones and Kang 1990; Kang and Jones 1995, 1997) 
models of CR transport. Kang and Jones (1995) using PPM hydrodynamics and Kang and Jones (1997) using TVD magnetohydro- 
dynamics methods, for example, provided comparisons between the CR transport approach applied here with our earlier "thermal 
leakage" injection model and Monte Carlo, hybrid plasma simulations, and in situ measurements of heliospheric shocks. While the 
newer injection scheme utilizes a significantly improved physical model to determine which particles are injected at a shock, both the 
new and the old underlying numerical schemes have much in common, as discussed in Gieseler et al. (2000). More recently Kang 
et al. (2001) provided convergence tests of the combined shock tracking, AMR, convection diffusion code upon which CRASH is 
based, again using our previous thermal leakage injection model. 

In order to extend the above test base to CRASH we present here several comparisons of solutions to results from these earlier 
works. In particular we have tested CRASH against piston-driven shock problems that were calculated with our PPM/CR code 
(Jones and Kang 1990; Kang, Jones and Ryu 1992). First, using a two-fluid version of our CRASH code with the modified entropy 
equation, we calculated a shock driven by a constant speed, plane piston with u p = 0.3 moving into a medium of p = 1, P g = 7 x 10 -4 , 
and P c = 3.5 x 10~ 4 . Closure parameters are assumed to be 7 C = 5/3 and < k >= 0.01. There is no injection of new CRs in this 
case, only reacceleration of the previous CR population. Thus, it is a test for nonlinear modified shocks of the dynamical coupling 
terms included in equations (2)-(3). Fig. 13 shows the evolution of the shock at t = 3, 6, and 9, which can be compared with Fig. 
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1 of Jones and Kang (1990). Agreement is excellent. The dashed horizontal line shows the steady state solution for P c , which is 
fully consistent with the CRASH simulation. This test calculation demonstrates that the two fluid version of the CRASH code can 
accurately follow the evolution of a strongly CR modified shock. The dynamical coupling between the CRs and the gas is the same 
in the two-fluid model as in the diffusion-convection, kinetic equation model, except for the multi-scale character introduced by a 
momentum dependent diffusion coefficient into the kinetic equation model. 

To extend the dynamical evaluation to include momentum dependent influences, we calculated another piston driven shock using 
the diffusion-convection version of CRASH. In this case the piston has velocity u p = l moving into a medium of p = 1, P g = P C = 0.001. 
Two simulations were done, one with and one without the modified entropy equation (5). The injection process was turned off and 
only the base-grid was used (i.e., / g max = 0). The diffusion coefficient was given the form n = p 5 . Fig. 14 shows the evolution of the 
shock at t = 10, 30, 50, 100 and 150, which was calculated with the modified entropy equation. The last of these compares directly 
with the dotted curves in Fig. 3 of Kang, Jones and Ryu (1992). Some details of the early evolution differ, which cause the exact 
locations and peaks of the transient density spike to differ somewhat, but the overall structures and the shock transitions at t = 150 are 
in very good agreement, including the compression, as well as the gas and CR pressures. The noted small early evolution differences 
among different codes have been seen before. They reflect the sensitivity to numerical details of the very rapid changes taking place 
in the shocks at the start of such simulations. 

Finally, we demonstrate the benefits of including the modified entropy equation in treatments of strongly modified CR shocks. 
Fig. 15 shows from two simulations the subshocks and precursor structure for the piston-driven shock shown in Fig. 14 at t = 20. 
The two sets of curves were calculated with the modified entropy condition enforced (S code: solid line) and without it (E code: 
dashed line). The lower left panel shows the specific entropy, \ogP g / p 1 . To the right of the subshock, where the flow is adiabatic, 
this quantity should be a constant. However, because the flow there is highly supersonic, the total energy is dominated by the kinetic 
energy. Consequently, errors in computing the gas pressures used in equations (2) and (3), as extracted from a combination of the 
total energy and momentum quantities lead to significant errors in the entropy of the preshock gas. That in turn, leads to errors in the 
gas subshock and CR acceleration behaviors. Enforcing entropy conservation in smooth flows eliminates the problem, as illustrated 
in the figure. 
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FIG. Al. — Thermal velocity, t>th, and "effective" injection velocity, «i n j, both normalized by the downstream flow velocity relative to the subshock, uj, are 
shown as functions of subshock Mach numbers. The effective injection velocity is also dependent on the inverse wave parameter, e. In the lower left panel, the 
ratios of Uj n j/uth are plotted for e = 0.2 — 0.35. Thermal leakage is faster for smaller v- m j/v t ^. In the lower right panel, the thermal distribution function, jm of the 
downstream gas for Mo = 100 shock is plotted against v/uj, where uj = u s /r sv fy, u s /c = 0.01, c is the speed of light, and is the subshock compression ratio. The 
transparency function, r esc , for the same shock is also shown for e = 0.3 (solid line) and 0.2 (dashed line). The axis on the left side is for log(gM), while the axis on 
the right side is for T esc . 
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FIG. A2. — Density profile in the refinement region around the shock for Mo = 100 and e = 0.2 model at t= 10 and 20 at the grid level l g = 1. 3, 5, and 7. 
are 200 zones at each grid level which are shown as filled points at f = 20 curves. 
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FIG. A3. — Time evolution of the Mq = 5 shock structure with e = 0.3 and /max = 5 refined grid levels. Gas density, pressure, flow velocity and CR pressure are 
shown at f/f = 0, 4, 8, 12, and 20. The right-facing shock drifts to the left, so the right most transitions correspond to the earliest time t = 0. 
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FIG. A5. — Time evolution of subshock velocity, u s , subshock Mach number, M s , preshock density, pi/po and postshock density, P2/P0 are shown for Mo = 2—30 
and e = 0.3. 





FIG. A6. — The ratio of total CR energy in the simulation box to the kinetic energy in the initial shock rest frame that has entered the simulation box from 
upstream, <3?(f), the postshock CR pressure in units of upstream ram pressure in the instantaneous shock frame, and time-averaged injection efficiency, £(f). Left 
three panels are for Mo = 3 - 30 and e = 0.3. Right three panels show the same quantities for Mo = 5 - 100 and e = 0.25. 
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FIG. A7. — Evolution of the CR distribution function at the shock, represented as g = p 4 f(p), is shown for the same shock as in Fig. 3 (Mo = 5 and e = 0.3). The 
transparency function is also plotted at t/t = and 20 for reference. The axis on the left side is for log(gM), while the axis on the right side is for r esc . 
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FIG. A8. — Evolution of the CR distribution function at the shock, represented as g = p 4 f(p), is shown for the same shock in Fig. 4 (Mo = 30 and e = 0.3). The 
transparency function is also plotted at t/t = and 20 for reference. The axis on the left side is for log(gM), while the axis on the right side is for r esc . 
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FIG. A9. — Gas density, pressure, CR adiabatic index and CR pressure are shown at t/t = 0, 10, 20, 50, 100, and 140 for a Mq = 200 shock. The maximum level 
of refinement was / max = 7. The shock drifts to the left, so the right most transitions correspond to the earliest time t = 0. 
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FIG. A10. — The ratio of total CR energy in the simulation box to the kinetic energy in the initial shock rest frame that has entered the simulation box from 
upstream, "5(f), the postshock CR pressure in units of upstream ram pressure in the instantaneous shock frame, and time-averaged injection efficiency, £(f). Left 
three panels are for Mq = 5-200 and e = 0.2. Right three panels are for Mq = 30 and e = 0.2, 0.23, and 0.25. 
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FIG. All. — Evolution of the CR distribution function at the shock, represented as g = p 4 f(p), is shown for the Mo = 100 shock with e = 0.2. The transparency 
function is also plotted at t/t a = and 100 for reference. The axis on the left side is for log(gM), while the axis on the right side is for T esc . 
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FIG. A 12. — Left three panels: preshock density, pi/ pa, postshock density, P2I pa, and subshock compression ratio for Mq = 5-200. For models with e = 0.3 and 
0.25 the numerical results at t = 20 are plotted, while for models with e = 0.2 the results at t = 100 are plotted. Right top panel: approximate time asymptotic values 
of <!> is plotted. The same symbols and line types are used as in pi/po plot. Right middle panel: time asymptotic values of the upstream flow velocity relative to the 
subshock, u' Q , at t = 100 are plotted for e = 0.2. Right bottom panel: time asymptotic postshock CR pressure normalized to the ram pressure of the upstream flow 

in the initial shock rest frame, P c .2/po«o ( sou d une ). an d to tne ram pressure of the upstream flow in the instantaneous shock rest frame, P c ,2/po«o 2 (dotted line), is 
plotted for e = 0.2. 
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FIG. A13. — A piston-driven shock with u p = 0.3, propagating into a uniform medium of p = l,P g = 7 X 10~ 4 ,P C = 3.5 X lff"*,7 c = 5/3, and < k >= 0.01. 
CR energy injection at the subshock is not included. The simulation was calculated with the two-fluid version of our CRASH code with the modified entropy (S) 
equation. The results are shown at t = 3, 6, and 9. The dashed line corresponds to the postshock P c in the steady state limit. 
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FIG. A 14. — A piston-driven shock with u p = 1., propagating into a uniform medium of p = l,P g =P C = 0.001, and ft = p° 5 . CR particle injection at the subshock 
was turned off. The simulation was calculated with the kinetic version of our CRASH code with the modified entropy (S) equation. The results are shown at 
r = 10,30,50, 100 and 150. 
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FIG. A15. — The same shock shown in Fig. 14 except at t = 20. The solid lines are for the simulation with the modified entropy equation, while the dashed lines 
are for the simulation with the total energy equation only. 



